#include <algorithm/analytical_model/algorithm/cpu/heliostat/RectangleHeliostat.h>
#include <algorithm/analytical_model/algorithm/cuda/GeometryUtils.cuh>
#include <algorithm/analytical_model/algorithm/cuda/Utils.cuh>
#include <algorithm/analytical_model/algorithm/cuda/VectorMath.cuh>
#include <algorithm/analytical_model/algorithm/cuda/heliostat/RectangleHeliostat.cuh>
#include <utils/Utils.h>

#include <bit>
#include <tuple>
#include <vector>

using namespace solar;
using namespace solar::cuda;
using namespace solar::cuda::inter_op;

namespace
{
    __global__ void mapMicrohelioCenterAndNormal(float3* d_microhelio_centers,
                                                 float3* d_microhelio_normals, float3 normal,
                                                 float3 helio_size, const int2 row_col,
                                                 const int2 sub_row_col, const float pixel_length,
                                                 const float2 gap, const float3 world_pos,
                                                 const size_t map_size)
    {
        int myId = getThreadID();
        // printf("线程　%llu\n",myId);
        if (myId >= map_size)
        {
            return;
        }

        int row = myId / (row_col.y * sub_row_col.y);
        int col = myId % (row_col.y * sub_row_col.y);

        int block_row = row / sub_row_col.x;
        int block_col = col / sub_row_col.y;

        // 1. Generate local micro-heliostats' centers
        d_microhelio_centers[myId].x =
            col * pixel_length + block_col * gap.x + pixel_length / 2 - helio_size.x / 2;
        d_microhelio_centers[myId].y = helio_size.y / 2;
        d_microhelio_centers[myId].z =
            row * pixel_length + block_row * gap.y + pixel_length / 2 - helio_size.z / 2;

        // 2. Generate micro-heliostats' normals
        d_microhelio_normals[myId] = normal;

        // 3. Transform local micro-heliostat center to world position.
        float3 local = d_microhelio_centers[myId];
        local = localToWorldRotate(local, normal);
        local = translate(local, world_pos);
        d_microhelio_centers[myId] = local;
        // printf("微定日镜row=%d col=%d center.x=%f y=%f z=%f 法向x=%f y=%f
        // z=%f\n",row,col,local.x,local.y,local.z,normal.x,normal.y,normal.z);
    }
} // namespace

void RectangleHeliostat::setNormal(const float3& focus_center, const float3& sunray_dir)
{
    this->focus_center_ = focus_center;
    // float3 local_center = ::make_float3(pos_.x, pos_.y, pos_.z);
    // float3 reflect_dir1 = focus_center - local_center;
    // float3 reflect_dir = normalize(reflect_dir1);
    // normal_ = normalize(reflect_dir - sunray_dir);
    structCopy(cpu::setNormal(structCopy<Float3>(pos_), structCopy<Float3>(focus_center),
                              structCopy<Float3>(sunray_dir)),
               normal_);
}

void RectangleHeliostat::setWorldVertex()
{
    // 左下 左上 右上 右下

    vertices_[0] = ::make_float3(-size_.x / 2, size_.y / 2, -size_.z / 2);
    vertices_[1] = ::make_float3(0, 0, size_.z) + vertices_[0];
    vertices_[2] = ::make_float3(size_.x, 0, size_.z) + vertices_[0];
    vertices_[3] = ::make_float3(size_.x, 0, 0) + vertices_[0];
    /*printf("局部定日镜　%f %f %f\n", vertices_[0].x, vertices_[0].y, vertices_[0].z);
    printf("局部定日镜　%f %f %f\n", vertices_[1].x, vertices_[1].y, vertices_[1].z);
    printf("局部定日镜　%f %f %f\n", vertices_[2].x, vertices_[2].y, vertices_[2].z);
    printf("局部定日镜　%f %f %f\n", vertices_[3].x, vertices_[3].y, vertices_[3].z);*/
    // First step: rotate
    vertices_[0] = localToWorldRotate(vertices_[0], normal_);
    vertices_[1] = localToWorldRotate(vertices_[1], normal_);
    vertices_[2] = localToWorldRotate(vertices_[2], normal_);
    vertices_[3] = localToWorldRotate(vertices_[3], normal_);
    // Second step: translate
    vertices_[0] = translate(vertices_[0], pos_);
    vertices_[1] = translate(vertices_[1], pos_);
    vertices_[2] = translate(vertices_[2], pos_);
    vertices_[3] = translate(vertices_[3], pos_);
    /*printf("定日镜normal　%f %f %f\n", normal_.x, normal_.y, normal_.z);
    printf("pos　%f %f %f\n", pos_.x, pos_.y, pos_.z);

    printf("定日镜　%f %f %f\n", vertices_[0].x, vertices_[0].y, vertices_[0].z);
    printf("定日镜　%f %f %f\n", vertices_[1].x, vertices_[1].y, vertices_[1].z);
    printf("定日镜　%f %f %f\n", vertices_[2].x, vertices_[2].y, vertices_[2].z);
    printf("定日镜　%f %f %f\n", vertices_[3].x, vertices_[3].y, vertices_[3].z);*/
}

void RectangleHeliostat::setNormalAndRotate(const float3& focus_center, const float3& sunray_dir)
{
    setNormal(focus_center, sunray_dir);
    setWorldVertex();
}

auto RectangleHeliostat::getDiscreteMicroHelioOriginsAndNormals(
    float3*& d_microhelio_centers, float3*& d_microhelio_normals) const -> int
{
    float2 subhelio_row_col_length;
    subhelio_row_col_length.x = (size_.z - gap_.y * (row_col_.x - 1)) / float(row_col_.x);
    subhelio_row_col_length.y = (size_.x - gap_.x * (row_col_.y - 1)) / float(row_col_.y);
    // cout << "subhelio_row_col_length----------" << subhelio_row_col_length.x << ", " <<
    // subhelio_row_col_length.y << endl;

    int2 sub_row_col;
    sub_row_col.x = int(subhelio_row_col_length.x / pixel_length_);
    sub_row_col.y = int(subhelio_row_col_length.y / pixel_length_);
    // cout << "sub_row_col------" << sub_row_col.x << ", " << sub_row_col.y << endl;

    int map_size = row_col_.x * row_col_.y * sub_row_col.x * sub_row_col.y;
    // cout << "map_size-----" << map_size << endl;

    int n_threads = 0;
    dim3 n_blocks;
    setThreadBlocks(n_blocks, n_threads, map_size);

    // Map micro-heliostat center  and normal in world position.
    if (d_microhelio_centers == nullptr)
    {
        CHECK_CUDA_ERRORS(cudaMalloc((void**)&d_microhelio_centers, sizeof(float3) * map_size));
    }
    if (d_microhelio_normals == nullptr)
    {
        CHECK_CUDA_ERRORS(cudaMalloc((void**)&d_microhelio_normals, sizeof(float3) * map_size));
    }

    mapMicrohelioCenterAndNormal<<<n_blocks, n_threads>>>(
        d_microhelio_centers, d_microhelio_normals, normal_, size_, row_col_, sub_row_col,
        pixel_length_, gap_, pos_, map_size);
    // printf("map_size: ");
    // printf("%d\t", map_size);
    return map_size;
}

void RectangleHeliostat::getSubHeliostatVertices(std::vector<float3>& subHeliostatVertexes) const
{
    subHeliostatVertexes.emplace_back(vertices_[0]);
    subHeliostatVertexes.emplace_back(vertices_[1]);
    subHeliostatVertexes.emplace_back(vertices_[2]);
}

auto RectangleHeliostat::getImagePlaneAxis() const -> std::tuple<float3, float3, float3>
{
#ifdef __CUDA_ARCH__
    float3 x_axis;
    float3 y_axis;
    float3 z_axis;
    float3 reflect_dir = focus_center_ - pos_;
    z_axis = normalize(-reflect_dir);
    float3 interv[4];
    #pragma unroll 4
    for (int i = 0; i < 4; i++)
    {
        std::tie(interv[i], std::ignore) =
            intersect(z_axis, focus_center_, vertices_[i], reflect_dir);
    }
    x_axis = normalize(interv[2] - interv[1]);
    y_axis = cross(z_axis, x_axis);
    return {x_axis, y_axis, z_axis};
#else
    const auto* ptr = std::bit_cast<const Float3*>(std::addressof(this->vertices_));
    auto result = solar::cpu::getImagePlaneAxis(structCopy<Float3>(focus_center_),
                                                structCopy<Float3>(pos_), ptr);
    auto& x_axis = std::get<0>(result);
    auto& y_axis = std::get<1>(result);
    auto& z_axis = std::get<2>(result);
    return {structCopy<float3>(x_axis), structCopy<float3>(y_axis), structCopy<float3>(z_axis)};
#endif
}
